##Study_2_TESS_analysis_and_figures.R
##analyses in paper and appendix
##plots marginal effects of regression tables
##R version 3.6.1 

##set working directory where files are located
setwd("....") ##SET WORKING DIRECTORY

library(foreign)

data1<- read.dta("Study_2_TESS_for_analysis.dta", convert.factors=FALSE)

##note that between setup and analysis, saving file converted . to _ in names
##analysis uses _

####################################
##            Analysis            ##
####################################
##Analyses drop 20 respondents who did not answer main 3 DVs (bill, congress, fairness)
##Note that in saving file in data set up, . was converted to _. Change back
summary(as.factor(data1$answer_main))

##Do people think the news covers this? (subsetting as with analyses) - referenced in conclusion
table(data1$Q8[data1$answer_main==1]) ##
(62+203)/(62+203+583+577+637) ##13% this nearly all or most
(583)/(62+203+583+577+637) ##28% this some of the time
(577+637)/(62+203+583+577+637) ##59% think few or none

#######################################
##  Sample demographics for appendix ##
#######################################
prop.table(table(data1$GENDER[data1$answer_main==1]))
prop.table(table(data1$AGE7[data1$answer_main==1])) ##top 2 categories combined
prop.table(table(data1$REGION4[data1$answer_main==1]))
prop.table(table(data1$RACETHNICITY[data1$answer_main==1])) ##categories 3 and 5 combined for other
prop.table(table(data1$pid3_indpluslean[data1$answer_main==1]))

######################################################
##  Balance checks (refernced but not in table      ##
######################################################
##Check for balance of demographics across conditions
##Balance using chi-squared
data.finish<- subset(data1, data1$answer_main==1)
library(MASS) ##(only run for this, causes issues with tdyr/dplyr, need to detach if)
##gender (also in SSI)
balance.gender<- table(data.finish$GENDER, data.finish$treatment_category)
prop.table(balance.gender, 2)
chisq.test(balance.gender)
##age (also in SSI)
balance.age<- table(data.finish$AGE7, data.finish$treatment_category)
prop.table(balance.age, 2)
chisq.test(balance.age)
##region (also in SSI)
balance.region<- table(data.finish$REGION4, data.finish$treatment_category)
prop.table(balance.region, 2)
chisq.test(balance.region)
##Race (also in SSI)
balance.race<- table(data.finish$RACETHNICITY, data.finish$treatment_category)
prop.table(balance.race, 2)
chisq.test(balance.race)
##Party 
balance.party<- table(data.finish$pid3_indpluslean, data.finish$treatment_category)
prop.table(balance.party, 2)
chisq.test(balance.party)
##education
balance.educ<- table(data.finish$EDUC4, data.finish$treatment_category)
prop.table(balance.educ, 2)
chisq.test(balance.educ)
##income
balance.income<- table(data.finish$INCOME, data.finish$treatment_category)
prop.table(balance.income, 2)
chisq.test(balance.income) 

detach("package:MASS", unload=TRUE)

#####################################
##    Time on Vignettes            ##
#####################################
##Control
summary(data1$duration_VIGNETTE[data1$answer_main==1 & data1$treatment_category=="Control"])
##Treatment
summary(data1$duration_VIGNETTE[data1$answer_main==1 & data1$treatment_category!="Control"])

######################################
##  Factual manipulation check      ##
######################################
table(data1$manip_correct[data1$answer_main==1]) ##1380 correct, 687 incorrect
1380/(1380+687)
######################################
##  Group level check of FMC        ##
######################################
##manipulation check correct
##Q8A = 1 = Rep,  2 = Dem, 3= Not specified

##if assigned to Democratic majority, what pct said majority was Democrats?
xtabs(~data1$maj_dem[data1$answer_main==1] + data1$Q8A[data1$answer_main==1])
372/(207+372+169) ##49.7% (main error was saying Republicans, which was ACTUAL current majority)
##if assigned to Republican majority, what pct said majority was Democrats?
xtabs(~data1$maj_rep[data1$answer_main==1] + data1$Q8A[data1$answer_main==1])
88/(1008+88+223) ##6.67%
##NOTE more in GOP majority because reandomized across all 8 conditions
##how much more likely to say Dems if in Dem majority condition?
0.4973262/0.06671721

######################################
##  Regressions                     ##
######################################

##Pool all respondents (including ind) - For Table 4
##congress
reg1.pool.cong<- lm(cong ~ treatment_bipart + treatment_min,
                    data=data1,
                    subset=c(data1$answer_main==1))
summary(reg1.pool.cong) 

##bill evaluation
reg1.pool.bill<- lm(bill ~ treatment_bipart + treatment_min,
                    data=data1,
                    subset=c(data1$answer_main==1))
summary(reg1.pool.bill) 

##feelings toward majority
reg1.pool.maj<- lm(ft_maj ~ treatment_bipart + treatment_min,
                   data=data1,
                   subset=c(data1$answer_main==1)) 
summary(reg1.pool.maj) 

##fairness of process
reg1.pool.fair<- lm(fair ~ treatment_bipart + treatment_min,
                    data=data1,
                    subset=c(data1$answer_main==1))
summary(reg1.pool.fair)


##calcs for substantive magnitude:
##bill
sd(data1$bill[data1$answer_main==1]) ##0.3049
##for bill eval, effect of bipartisan treatment is -0.126
0.126/0.3049 ##0.413
##minority is -0.13451
-0.13451/0.3049 ##-0.44

##cong
sd(data1$cong[data1$answer_main==1]) 
##bipart
-0.03256/0.1984859 ##-0.16
##minority
-0.02660/0.1984859 ##-0.13

##fair
sd(data1$fair[data1$answer_main==1]) 
##bipart
-0.15018/0.2783247 ##0.54
##minority
-0.15254/0.2783247 ##0.55

##ft_maj
sd(data1$ft_maj[data1$answer_main==1], na.rm=TRUE) 
##bipart
-0.04658/0.313137 ##-0.15
##minority
-0.05154/0.313137 ##-0.16

###############################
##Interaction of treatments with respondents who align with the opposing party, exclude independents - Table 5

##bill evaluation
reg3.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
               data=data1,
               subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.bill) 

##congress
reg3.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
               data=data1,
               subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.cong) 

##fairness of process
reg3.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
               data=data1,
               subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.fair)

##feelings toward majority
reg3.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
              data=data1,
              subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.maj) 


##substantive magnitude among partisans
##bill
sd(data1$bill[data1$answer_main==1 & data1$pid3!="Ind"])
##bipart
-0.11966/0.3143588 ##-0.38
##minority
-0.09389/0.3143588 ##-0.30

##fair
sd(data1$fair[data1$answer_main==1 & data1$pid3!="Ind"])
##bipart
-0.13798/0.2831842 ##-0.487
##minority
-0.09219/0.2831842##-0.3255

############################
##motivations, interactions with minority voters, exlcudes pure independents -Table 6
reg3.motive.defeat<- lm(motive_defeat ~ treatment_bipart*maj_other + treatment_min*maj_other,
                        data=data1,
                        subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.defeat) 

reg3.motive.win<- lm(motive_win ~ treatment_bipart*maj_other + treatment_min*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.win) 

reg3.motive.refuselisten<- lm(motive_refuselisten ~ treatment_bipart*maj_other + treatment_min*maj_other,
                              data=data1,
                              subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.refuselisten) 

reg3.motive.nofair<- lm(motive_nofair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                        data=data1,
                        subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.nofair) 

reg3.motive.goodpolicy<- lm(motive_goodpolicy ~ treatment_bipart*maj_other + treatment_min*maj_other,
                            data=data1,
                            subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.goodpolicy) 

reg3.motive.base<- lm(motive_base ~ treatment_bipart*maj_other + treatment_min*maj_other,
                      data=data1,
                      subset=c(data1$pid3!="Ind" & data1$answer_main==1))
summary(reg3.motive.base)  


#########################################
##    Appendix Tables                  ##
#########################################
##Robustness checks for Table 5 for appendix 

##manipulation check correct Table E5
##Pool all conditions, exclude independents
##bill evaluation
reg3.manip.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$answer_main==1 & manip_correct==1))
summary(reg3.manip.bill) 

##congress
reg3.manip.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$answer_main==1 & manip_correct==1))
summary(reg3.manip.cong) 

##feelings toward majority
reg3.manip.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
                    data=data1,
                    subset=c(data1$pid3!="Ind" & data1$answer_main==1 & manip_correct==1))
summary(reg3.manip.maj) 

##fairness of process
reg3.manip.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$answer_main==1 & manip_correct==1))
summary(reg3.manip.fair)


##############################
##subset to strong partisans Table E6
##bill evaluation
reg3.strongpart.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
                          data=data1,
                          subset=c(data1$pid3!="Ind" & data1$answer_main==1 & strong_partisan==1))
summary(reg3.strongpart.bill) 

##congress
reg3.strongpart.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
                          data=data1,
                          subset=c(data1$pid3!="Ind" & data1$answer_main==1 & strong_partisan==1))
summary(reg3.strongpart.cong) 

##feelings toward majority
reg3.strongpart.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
                         data=data1,
                         subset=c(data1$pid3!="Ind" & data1$answer_main==1 & strong_partisan==1))
summary(reg3.strongpart.maj) 

##fairness of process
reg3.strongpart.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                          data=data1,
                          subset=c(data1$pid3!="Ind" & data1$answer_main==1 & strong_partisan==1))
summary(reg3.strongpart.fair)


#################################
##subset to high knowledge (median or more, 2+) Table E7
##bill evaluation
reg3.highknow.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
                        data=data1,
                        subset=c(data1$pid3!="Ind" & data1$answer_main==1 & knowledge_correct>=2))
summary(reg3.highknow.bill) 

##congress
reg3.highknow.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
                        data=data1,
                        subset=c(data1$pid3!="Ind" & data1$answer_main==1 & knowledge_correct>=2))
summary(reg3.highknow.cong) 

##feelings toward majority
reg3.highknow.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
                       data=data1,
                       subset=c(data1$pid3!="Ind" & data1$answer_main==1 & knowledge_correct>=2))
summary(reg3.highknow.maj) 

##fairness of process
reg3.highknow.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                        data=data1,
                        subset=c(data1$pid3!="Ind" & data1$answer_main==1 & knowledge_correct>=2))
summary(reg3.highknow.fair)


#################################
##Just independents (pure and leaning) - Appendix Table E8
##bill evaluation
reg1.ind.bill<- lm(bill ~ treatment_bipart + treatment_min,
                   data=data1,
                   subset=c(data1$pid3_indpluslean=="Ind" & data1$answer_main==1))
summary(reg1.ind.bill) 

##congress
reg1.ind.cong<- lm(cong ~ treatment_bipart + treatment_min,
                   data=data1,
                   subset=c(data1$pid3_indpluslean=="Ind" & data1$answer_main==1))
summary(reg1.ind.cong) 

##feelings toward majority
reg1.ind.maj<- lm(ft_maj ~ treatment_bipart + treatment_min,
                  data=data1,
                  subset=c(data1$pid3_indpluslean=="Ind" & data1$answer_main==1))
summary(reg1.ind.maj) 

##fairness of process
reg1.ind.fair<- lm(fair ~ treatment_bipart + treatment_min,
                   data=data1,
                   subset=c(data1$pid3_indpluslean=="Ind" & data1$answer_main==1))
summary(reg1.ind.fair)


################################
##Split party/issue conditions - Appendix Table E9

##bill green energy alternative
reg2.bill.green<- lm(bill ~ treatment_bipart_green*maj_other + treatment_min_green*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax==0 & data1$answer_main==1))
summary(reg2.bill.green)

##congress green energy alternative
reg2.cong.green<- lm(cong ~ treatment_bipart_green*maj_other + treatment_min_green*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax==0 & data1$answer_main==1))
summary(reg2.cong.green)

##fairness of process
reg2.fair.green<- lm(fair ~ treatment_bipart_green*maj_other + treatment_min_green*maj_other,
                     data=data1,
                     subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax==0 & data1$answer_main==1))
summary(reg2.fair.green)

############
##bill gastax alternative
reg2.bill.gastax<- lm(bill ~ treatment_bipart_gastax*maj_other + treatment_min_gastax*maj_other,
                      data=data1,
                      subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax_plus_control==1 & data1$answer_main==1))
summary(reg2.bill.gastax) 

##congress gas tax alternative
reg2.cong.gastax<- lm(cong ~ treatment_bipart_gastax*maj_other + treatment_min_gastax*maj_other,
                      data=data1,
                      subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax_plus_control==1 & data1$answer_main==1))
summary(reg2.cong.gastax) 

##fairness of process
reg2.fair.gastax<- lm(fair ~ treatment_bipart_gastax*maj_other + treatment_min_gastax*maj_other,
                      data=data1,
                      subset=c(data1$pid3!="Ind" & data1$maj_rep==1 & data1$gastax_plus_control==1 & data1$answer_main==1))
summary(reg2.fair.gastax)

#############
## Dem majority (alternative oil)
##bill 
reg2.bill.oil<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
                   data=data1,
                   subset=c(data1$pid3!="Ind" & data1$maj_dem==1 & data1$answer_main==1))
summary(reg2.bill.oil) 

##congress 
reg2.cong.oil<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
                   data=data1,
                   subset=c(data1$pid3!="Ind" & data1$maj_dem==1 & data1$answer_main==1))
summary(reg2.cong.oil) 

##fairness of process
reg2.fair.oil<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                   data=data1,
                   subset=c(data1$pid3!="Ind" & data1$maj_dem==1 & data1$answer_main==1))
summary(reg2.fair.oil)

#############################
##Interaction with gas tax Appendix footnote 11 and Appendix Table E10
##Pool all Republican majority conditions, exclude independents - add 3-way interaction with gas tax
##bill evaluation
reg3.gasinteract.bill<- lm(bill ~ treatment_bipart*maj_other*gastax + treatment_min*maj_other*gastax -(gastax) - (maj_other:gastax),
                           data=data1,
                           subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg3.gasinteract.bill) 

##congress
reg3.gasinteract.cong<- lm(cong ~ treatment_bipart*maj_other*gastax + treatment_min*maj_other*gastax -(gastax) - (maj_other:gastax),
                           data=data1,
                           subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg3.gasinteract.cong) 

##feelings toward majority
reg3.gasinteract.maj<- lm(ft_maj ~ treatment_bipart*maj_other*gastax + treatment_min*maj_other*gastax -(gastax) - (maj_other:gastax),
                          data=data1,
                          subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg3.gasinteract.maj) 

##fairness of process
reg3.gasinteract.fair<- lm(fair ~ treatment_bipart*maj_other*gastax + treatment_min*maj_other*gastax -(gastax) - (maj_other:gastax),
                           data=data1,
                           subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg3.gasinteract.fair)


#######################################################
##    Robustness check on Republican majority only   ##
#######################################################
##Because most common error in FMC was saying Republican majorty, 
##look at only Republican majority treatments where much higher fraction got manipulation check correct
xtabs(~data1$maj_rep + data1$Q8A)
1014/(1014+91+225)

###############################
##Appendix footnote 8 check among only respondents assigned to Republican majority
##bill evaluation
reg.gop.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
                  data=data1,
                  subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg.gop.bill) 

##congress
reg.gop.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
                  data=data1,
                  subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg.gop.cong) 

##feelings toward majority
reg.gop.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
                 data=data1,
                 subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg.gop.maj) 

##fairness of process
reg.gop.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
                  data=data1,
                  subset=c(data1$pid3!="Ind" & data1$answer_main==1 & data1$maj_rep==1))
summary(reg.gop.fair) 



#######################################################
##           Output regressions                      ##
#######################################################
source("regtable.R")

##set working directory for output
setwd(".../output") ##SET WORKING DIRECTORY HERE

##Pooled models, no interaction - Table 4
outtable.rtf(list("(1) Bill"=reg1.pool.bill, "(2) Congress"=reg1.pool.cong, "(3) Majority"=reg1.pool.maj, "(4) Fair Process"=reg1.pool.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Table 4.rtf")

##Models with all conditions, excluding independents - Table 5
outtable.rtf(list("(1) Bill"=reg3.bill, "(2) Congress"=reg3.cong, "(3) Majority"=reg3.maj, "(4) Fair Process"=reg3.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Table 5.rtf")

##Motives Models with all conditions, excluding independents - Table 6
outtable.rtf(list("(1) Defeat"=reg3.motive.defeat, "(2) Win"=reg3.motive.win,   
                  "(3) Refuse Listen"=reg3.motive.refuselisten, "(4) No Fair"=reg3.motive.nofair,
                  "(4) Good Policy"=reg3.motive.goodpolicy, "(6) Base"=reg3.motive.base),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
            "Table 6.rtf")

##Just independents - Appendix Table E8
outtable.rtf(list("(1) Bill"=reg1.ind.bill, "(2) Congress"=reg1.ind.cong, "(3) Majority"=reg1.ind.maj, "(4) Fair Process"=reg1.ind.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E8.rtf")


##Models without Ind, broken by majority/issue  - Appendix E9
outtable.rtf(list("(1) Bill (GOP/Green)"=reg2.bill.green, "(2) Bill (GOP/Gas)"=reg2.bill.gastax, "(3) Bill (Dem/Oil)"=reg2.bill.oil,
                  "(4) Congress (GOP/Green)"=reg2.cong.green, "(5) Congress (GOP/Gas)"=reg2.cong.gastax, "(6) Congress (Dem/Oil)"=reg2.cong.oil,
                  "(7) Fair Process (GOP/Green)"=reg2.fair.green, "(8) Fair Process (GOP/Gas)"=reg2.fair.gastax, "(9) Fair Process (Dem/Oil)"=reg2.fair.oil),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters"),
                              c("treatment_bipart_green", "Ignore Bipartisan"),
                              c("treatment_min_green", "Ignore Minority"),
                              c("treatment_bipart_green:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min_green", "Minority X Minority Voters"),
                              c("treatment_bipart_gastax", "Ignore Bipartisan"),
                              c("treatment_min_gastax", "Ignore Minority"),
                              c("treatment_bipart_gastax:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min_gastax", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E9.rtf")


##Appendix manipulation check correct - Appendix Table E5
outtable.rtf(list("(1) Bill"=reg3.manip.bill, "(2) Congress"=reg3.manip.cong, "(3) Majority"=reg3.manip.maj, "(4) Fair Process"=reg3.manip.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E5.rtf")

##Appendix strong partisans - Appendix Table E6
outtable.rtf(list("(1) Bill"=reg3.strongpart.bill, "(2) Congress"=reg3.strongpart.cong, "(3) Majority"=reg3.strongpart.maj, "(4) Fair Process"=reg3.strongpart.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E6.rtf")

##Appendix high knowledge -Appendix TAble E7
outtable.rtf(list("(1) Bill"=reg3.highknow.bill, "(2) Congress"=reg3.highknow.cong, "(3) Majority"=reg3.highknow.maj, "(4) Fair Process"=reg3.highknow.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E7.rtf")

##Appendix gas tax interaction -Models with Republican majority conditions, excluding independents - appendix Table E10
##Note that manually exlcudes lower order interactions
outtable.rtf(list("(1) Bill"=reg3.gasinteract.bill, "(2) Congress"=reg3.gasinteract.cong, "(3) Majority"=reg3.gasinteract.maj, "(4) Fair Process"=reg3.gasinteract.fair),
             replacelist=list(c("(Intercept)", "Constant"),
                              c("treatment_bipart", "Ignore Bipartisan"),
                              c("treatment_min", "Ignore Minority"),
                              c("maj_other", "Minority Voters"),
                              c("treatment_bipart:maj_other", "Bipartisan X Minority Voters"),
                              c("maj_other:treatment_min", "Minority X Minority Voters"),
                              c("treatment_bipart:gastax", "Bipartisan X Gas Tax"),
                              c("gastax:treatment_min", "Minority X Gas Tax"),
                              c("treatment_bipart:maj_other:gastax", "Bipartisan X Min. Voters X Gas Tax"),
                              c("maj_other:gastax:treatment_min", "Minority X Min. Voters X Gas Tax")),
             p.levels =c(0.10,0.05,0.01,0.001),
             scientific = 5,
             digits = 3,
             p.levels.labels=c("^", "*","**","***"),
             "Appendix Table E10.rtf")


################################################
##Plot marginal effects                       ##
################################################

##Analysis with dplyr
library(dplyr)
library(tidyr)
library(infer)
library(ggplot2)
library(gridExtra)
library(ggeffects)
library(margins)

##marginal effect plot resources
#https://cran.r-project.org/web/packages/sjPlot/vignettes/plot_marginal_effects.html
#https://strengejacke.github.io/ggeffects/
####For a resource, see: http://www.statsblogs.com/2013/08/27/creating-marginal-effect-plots-for-linear-regression-models-in-r/, which is the R version of the Brambor, Golder, et al. resource: http://mattgolder.com/files/interactions/interaction1.pdf

###############################
##Pool all conditions, exclude independents - Table 4
##bill evaluation
summary(reg3.bill) 

##congress
summary(reg3.cong) 

##feelings toward majority
summary(reg3.maj) 

##fairness of process
summary(reg3.fair)

#############################################
##Marginal effects
##Manual and use tibble to get to ggplot
##Bill
m.bill<- lm(bill ~ treatment_bipart*maj_other + treatment_min*maj_other,
            data=data1,
            subset=c(data1$pid3!="Ind" & data1$answer_main==1))

# pull out coefficient estimates
bill.beta.hat <- coef(m.bill)
bill.beta.hat

# pull out the covariance matrix
bill.cov <- vcov(m.bill)
bill.cov

# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(data1$maj_other), max(data1$maj_other), length.out = 2)
z0

##1) ignore bipartisan
# calculate the instantaneous effect of x as z varies
bill.dy.dx.bipart <- bill.beta.hat["treatment_bipart"] + bill.beta.hat["treatment_bipart:maj_other"]*z0
bill.dy.dx.bipart

# calculate the standard error of each estimated effect
bill.se.dy.dx.bipart <- sqrt(bill.cov["treatment_bipart", "treatment_bipart"] + 
                   z0^2*bill.cov["treatment_bipart:maj_other", "treatment_bipart:maj_other"] + 
                   2*z0*bill.cov["treatment_bipart", "treatment_bipart:maj_other"])
bill.se.dy.dx.bipart

# calculate upper and lower bounds of a 95% CI 
bill.upr.bipart <- bill.dy.dx.bipart + 1.96*bill.se.dy.dx.bipart
bill.lwr.bipart <- bill.dy.dx.bipart - 1.96*bill.se.dy.dx.bipart

##2) ignore minority
# calculate the instantaneous effect of x as z varies
bill.dy.dx.min <- bill.beta.hat["treatment_min"] + bill.beta.hat["maj_other:treatment_min"]*z0
bill.dy.dx.min

# calculate the standard error of each estimated effect
bill.se.dy.dx.min <- sqrt(bill.cov["treatment_min", "treatment_min"] + 
                          z0^2*bill.cov["maj_other:treatment_min", "maj_other:treatment_min"] + 
                          2*z0*bill.cov["treatment_min", "maj_other:treatment_min"])
bill.se.dy.dx.min

# calculate upper and lower bounds of a 95% CI 
bill.upr.min <- bill.dy.dx.min + 1.96*bill.se.dy.dx.min
bill.lwr.min <- bill.dy.dx.min - 1.96*bill.se.dy.dx.min

##combine as tibble
margins.bill<- tibble(treatment=c("Ignore Bipartisan", "Ignore Bipartisan", "Ignore Minority", "Ignore Minority"), 
                                 party=c("Majority Voters", "Minority Voters", "Majority Voters", "Minority Voters"), 
                                 y=c(bill.dy.dx.bipart, bill.dy.dx.min), se=c(bill.se.dy.dx.bipart, bill.se.dy.dx.min), 
                                 upr=c(bill.upr.bipart, bill.upr.min), lwr=c(bill.lwr.bipart, bill.lwr.min))
margins.bill

fig3a<- ggplot(margins.bill) +
  geom_errorbar(aes(x=party, ymin=lwr, ymax=upr), ##95% CI
                width=.15) +
  geom_point(aes(x=party, y=y), size=2) + # add the points on
  facet_grid(~treatment) + # split into plots, side by side
  scale_y_continuous(limits=c(-.25,.05)) + 
  theme(axis.text=element_text(size=10), axis.title.y = element_text(size = 10),
        plot.title = element_text(size=10))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ggtitle("Dependent Variable: Bill Support") + ylab("Marginal Effect of Treatment") +
  geom_hline(yintercept=0, color="grey")


##############################################
##cong
m.cong<- lm(cong ~ treatment_bipart*maj_other + treatment_min*maj_other,
            data=data1,
            subset=c(data1$pid3!="Ind" & data1$answer_main==1))

# pull out coefficient estimates
cong.beta.hat <- coef(m.cong)
cong.beta.hat

# pull out the covariance matrix
cong.cov <- vcov(m.cong)
cong.cov

# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(data1$maj_other), max(data1$maj_other), length.out = 2)
z0

##1) ignore bipartisan
# calculate the instantaneous effect of x as z varies
cong.dy.dx.bipart <- cong.beta.hat["treatment_bipart"] + cong.beta.hat["treatment_bipart:maj_other"]*z0
cong.dy.dx.bipart

# calculate the standard error of each estimated effect
cong.se.dy.dx.bipart <- sqrt(cong.cov["treatment_bipart", "treatment_bipart"] + 
                               z0^2*cong.cov["treatment_bipart:maj_other", "treatment_bipart:maj_other"] + 
                               2*z0*cong.cov["treatment_bipart", "treatment_bipart:maj_other"])
cong.se.dy.dx.bipart

# calculate upper and lower bounds of a 95% CI 
cong.upr.bipart <- cong.dy.dx.bipart + 1.96*cong.se.dy.dx.bipart
cong.lwr.bipart <- cong.dy.dx.bipart - 1.96*cong.se.dy.dx.bipart

##2) ignore minority
# calculate the instantaneous effect of x as z varies
cong.dy.dx.min <- cong.beta.hat["treatment_min"] + cong.beta.hat["maj_other:treatment_min"]*z0
cong.dy.dx.min

# calculate the standard error of each estimated effect
cong.se.dy.dx.min <- sqrt(cong.cov["treatment_min", "treatment_min"] + 
                            z0^2*cong.cov["maj_other:treatment_min", "maj_other:treatment_min"] + 
                            2*z0*cong.cov["treatment_min", "maj_other:treatment_min"])
cong.se.dy.dx.min

# calculate upper and lower bounds of a 95% CI 
cong.upr.min <- cong.dy.dx.min + 1.96*cong.se.dy.dx.min
cong.lwr.min <- cong.dy.dx.min - 1.96*cong.se.dy.dx.min

##combine as tibble
margins.cong<- tibble(treatment=c("Ignore Bipartisan", "Ignore Bipartisan", "Ignore Minority", "Ignore Minority"), 
                      party=c("Majority Voters", "Minority Voters", "Majority Voters", "Minority Voters"), 
                      y=c(cong.dy.dx.bipart, cong.dy.dx.min), se=c(cong.se.dy.dx.bipart, cong.se.dy.dx.min), 
                      upr=c(cong.upr.bipart, cong.upr.min), lwr=c(cong.lwr.bipart, cong.lwr.min))
margins.cong

fig3b<- ggplot(margins.cong) +
  geom_errorbar(aes(x=party, ymin=lwr, ymax=upr), ##95% CI
                width=.15) +
  geom_point(aes(x=party, y=y), size=2) + # add the points on
  facet_grid(~treatment) + # split into plots, side by side
  scale_y_continuous(limits=c(-.25,.05)) + 
  theme(axis.text=element_text(size=10), axis.title.y = element_text(size = 10),
        plot.title = element_text(size=10))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ggtitle("Dependent Variable: Confidence in Congress") + ylab("Marginal Effect of Treatment") +
  geom_hline(yintercept=0, color="grey")


##############################################
##maj
m.maj<- lm(ft_maj ~ treatment_bipart*maj_other + treatment_min*maj_other,
            data=data1,
            subset=c(data1$pid3!="Ind" & data1$answer_main==1))

# pull out coefficient estimates
maj.beta.hat <- coef(m.maj)
maj.beta.hat

# pull out the covariance matrix
maj.cov <- vcov(m.maj)
maj.cov

# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(data1$maj_other), max(data1$maj_other), length.out = 2)
z0

##1) ignore bipartisan
# calculate the instantaneous effect of x as z varies
maj.dy.dx.bipart <- maj.beta.hat["treatment_bipart"] + maj.beta.hat["treatment_bipart:maj_other"]*z0
maj.dy.dx.bipart

# calculate the standard error of each estimated effect
maj.se.dy.dx.bipart <- sqrt(maj.cov["treatment_bipart", "treatment_bipart"] + 
                               z0^2*maj.cov["treatment_bipart:maj_other", "treatment_bipart:maj_other"] + 
                               2*z0*maj.cov["treatment_bipart", "treatment_bipart:maj_other"])
maj.se.dy.dx.bipart

# calculate upper and lower bounds of a 95% CI 
maj.upr.bipart <- maj.dy.dx.bipart + 1.96*maj.se.dy.dx.bipart
maj.lwr.bipart <- maj.dy.dx.bipart - 1.96*maj.se.dy.dx.bipart

##2) ignore minority
# calculate the instantaneous effect of x as z varies
maj.dy.dx.min <- maj.beta.hat["treatment_min"] + maj.beta.hat["maj_other:treatment_min"]*z0
maj.dy.dx.min

# calculate the standard error of each estimated effect
maj.se.dy.dx.min <- sqrt(maj.cov["treatment_min", "treatment_min"] + 
                            z0^2*maj.cov["maj_other:treatment_min", "maj_other:treatment_min"] + 
                            2*z0*maj.cov["treatment_min", "maj_other:treatment_min"])
maj.se.dy.dx.min

# calculate upper and lower bounds of a 95% CI 
maj.upr.min <- maj.dy.dx.min + 1.96*maj.se.dy.dx.min
maj.lwr.min <- maj.dy.dx.min - 1.96*maj.se.dy.dx.min

##combine as tibble
margins.maj<- tibble(treatment=c("Ignore Bipartisan", "Ignore Bipartisan", "Ignore Minority", "Ignore Minority"), 
                      party=c("Majority Voters", "Minority Voters", "Majority Voters", "Minority Voters"), 
                      y=c(maj.dy.dx.bipart, maj.dy.dx.min), se=c(maj.se.dy.dx.bipart, maj.se.dy.dx.min), 
                      upr=c(maj.upr.bipart, maj.upr.min), lwr=c(maj.lwr.bipart, maj.lwr.min))
margins.maj

fig3c<- ggplot(margins.maj) +
  geom_errorbar(aes(x=party, ymin=lwr, ymax=upr), ##95% CI
                width=.15) +
  geom_point(aes(x=party, y=y), size=2) + # add the points on
  facet_grid(~treatment) + # split into plots, side by side
  scale_y_continuous(limits=c(-.25,.05)) + 
  theme(axis.text=element_text(size=10), axis.title.y = element_text(size = 10),
        plot.title = element_text(size=10))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ggtitle("Dependent Variable: Feelings Toward Majority Party") + ylab("Marginal Effect of Treatment") +
  geom_hline(yintercept=0, color="grey")


##############################################
##fair
m.fair<- lm(fair ~ treatment_bipart*maj_other + treatment_min*maj_other,
            data=data1,
            subset=c(data1$pid3!="Ind" & data1$answer_main==1))

# pull out coefficient estimates
fair.beta.hat <- coef(m.fair)
fair.beta.hat

# pull out the covariance matrix
fair.cov <- vcov(m.fair)
fair.cov

# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(data1$maj_other), max(data1$maj_other), length.out = 2)
z0

##1) ignore bipartisan
# calculate the instantaneous effect of x as z varies
fair.dy.dx.bipart <- fair.beta.hat["treatment_bipart"] + fair.beta.hat["treatment_bipart:maj_other"]*z0
fair.dy.dx.bipart

# calculate the standard error of each estimated effect
fair.se.dy.dx.bipart <- sqrt(fair.cov["treatment_bipart", "treatment_bipart"] + 
                               z0^2*fair.cov["treatment_bipart:maj_other", "treatment_bipart:maj_other"] + 
                               2*z0*fair.cov["treatment_bipart", "treatment_bipart:maj_other"])
fair.se.dy.dx.bipart

# calculate upper and lower bounds of a 95% CI 
fair.upr.bipart <- fair.dy.dx.bipart + 1.96*fair.se.dy.dx.bipart
fair.lwr.bipart <- fair.dy.dx.bipart - 1.96*fair.se.dy.dx.bipart

##2) ignore minority
# calculate the instantaneous effect of x as z varies
fair.dy.dx.min <- fair.beta.hat["treatment_min"] + fair.beta.hat["maj_other:treatment_min"]*z0
fair.dy.dx.min

# calculate the standard error of each estimated effect
fair.se.dy.dx.min <- sqrt(fair.cov["treatment_min", "treatment_min"] + 
                            z0^2*fair.cov["maj_other:treatment_min", "maj_other:treatment_min"] + 
                            2*z0*fair.cov["treatment_min", "maj_other:treatment_min"])
fair.se.dy.dx.min

# calculate upper and lower bounds of a 95% CI 
fair.upr.min <- fair.dy.dx.min + 1.96*fair.se.dy.dx.min
fair.lwr.min <- fair.dy.dx.min - 1.96*fair.se.dy.dx.min

##combine as tibble
margins.fair<- tibble(treatment=c("Ignore Bipartisan", "Ignore Bipartisan", "Ignore Minority", "Ignore Minority"), 
                      party=c("Majority Voters", "Minority Voters", "Majority Voters", "Minority Voters"), 
                      y=c(fair.dy.dx.bipart, fair.dy.dx.min), se=c(fair.se.dy.dx.bipart, fair.se.dy.dx.min), 
                      upr=c(fair.upr.bipart, fair.upr.min), lwr=c(fair.lwr.bipart, fair.lwr.min))
margins.fair

fig3d<- ggplot(margins.fair) +
  geom_errorbar(aes(x=party, ymin=lwr, ymax=upr), ##95% CI
                width=.15) +
  geom_point(aes(x=party, y=y), size=2) + # add the points on
  facet_grid(~treatment) + # split into plots, side by side
  scale_y_continuous(limits=c(-.25,.05)) + 
  theme(axis.text=element_text(size=10), axis.title.y = element_text(size = 10),
        plot.title = element_text(size=10))+
  theme(axis.ticks.x=element_blank(), axis.title.x=element_blank(),
        #axis.title.y=element_blank(),
        strip.text.x = element_text(size=10)) + # how to change facet label font size
  ggtitle("Dependent Variable: Fairness of Process") + ylab("Marginal Effect of Treatment") +
  geom_hline(yintercept=0, color="grey")


#######################################
##combine plots
grid.arrange(fig3a, fig3b, fig3c, fig3d, nrow=2)
fig3.g<- arrangeGrob(fig3a, fig3b, fig3c, fig3d, nrow=2)
ggsave("figure 3_study 2 interactions.tiff", fig3.g, height=10, width=12, dpi=350)



